use micro-macro-2017, replace 

capture program drop vpc3
program vpc3, eclass 




local totvar = exp(_b[lns1_1_1:_cons])^2 + exp(_b[lns2_1_1:_cons])^2 + exp(_b[lnsig_e:_cons])^2
local districtvar = exp(_b[lns1_1_1:_cons])^2
local localvar = exp(_b[lns2_1_1:_cons])^2
local indivar =  exp(_b[lnsig_e:_cons])^2
display "% district: " `districtvar' / `totvar' * 100
display "% locality: " `localvar' / `totvar' * 100
display "% individual: " `indivar' / `totvar' * 100
matrix vpc = (0,0,0)
matrix vpc[1,1] = `indivar' / `totvar' * 100
matrix vpc[1,2] = `localvar' / `totvar' * 100
matrix vpc[1,3] = `districtvar' / `totvar' * 100
matrix colnames vpc = individual local county 
ereturn matrix vpc = vpc
end 


*********************
* Some empty models *
*********************

* Initialise matrix to collect results
matrix emptyvpcs = (0,0,0) \ (0,0,0) \ (0,0,0) \ (0,0,0) \ (0,0,0) \ (0,0,0) 
matrix colnames emptyvpcs = Individual Locality County
* This is annoying, because subsetting in Stata works only for extraction 

* Cultural threat
mixed cultthreat || district: || location:
est store emptycultthreat 
vpc3

matrix emptyvpcs = e(vpc)

* For completeness's sake
mixed econthreat || district: || location: 
est store emptyeconthreat
vpc3 

* Islamophobia
mixed islamophobia || district: || location: 
est store emptyislamophobia
vpc3 

matrix emptyvpcs = emptyvpcs \e(vpc)

* Populism. Difficult to get the empty model to converge 
mixed populism || district: || location: , matlog
est store emptypopulism
vpc3 

matrix emptyvpcs = emptyvpcs \e(vpc)

* RWA 
mixed rwa || district: || location: 
est store emptyrwa
vpc3 

matrix emptyvpcs = emptyvpcs \e(vpc)

* Conventionalism
mixed conventionalism  || district: || location: 
est store emptyconvention 
vpc3 

matrix emptyvpcs = emptyvpcs \e(vpc)

* Regional exclusion 
mixed regexclusion  || district: || location: 
est store emptyregexclusion  
vpc3 

matrix emptyvpcs = emptyvpcs \e(vpc)

* Localism 
mixed localism  || district: || location: 
est store emptylocalism   
vpc3 
matrix emptyvpcs = emptyvpcs \e(vpc)

* Set rownames 
matrix rownames emptyvpcs = "Cultural Threat" "Islamophobia" "Populism" "Auth. submission/aggression" "Conventionalism" "Regional exclusion" "Localism"

* Output
esttab matrix(emptyvpcs, fmt(1)) using emptyvpctable.rtf , replace mtitles("Variance")


*******************************
* Models with micro variables *
*******************************



* Ok. run a substantive model for cultural threat including gender, class, education etc.
mixed cultthreat i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion || district: || location: 
est store cultthreat_no_contextuals 

* Repeat for RWA
mixed rwa i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion || district: || location: 
est store rwa_no_contextuals 

* And for populism 
mixed populism i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion || district: || location: 
est store populism_no_contextuals 

* Islamophobia
mixed islamophobia i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion || district: || location: 
est store islamophobia_no_contextuals 

* Regional Exclusion 
mixed regexclusion i.male i.higheduc i.agecat i.worker i.unemployedsick  || district: || location: 
est store regexclusion_no_contextuals 

* Localism 
mixed localism i.male i.higheduc i.agecat i.worker i.unemployedsick  || district: || location: 
est store localism_no_contextuals 


* Ok. How much does including the individual-level variables reduce contextual variance?

foreach model in cultthreat islamophobia populism rwa regexclusion localism {
        display "Dependent variable: `model'"

        qui est restore empty`model'
        local emptycounty = exp(_b[lns1_1_1:_cons])^2
        local emptylocal = exp(_b[lns2_1_1:_cons])^2
        qui est restore `model'_no_contextuals 

        local individualcounty = exp(_b[lns1_1_1:_cons])^2
        local individuallocal = exp(_b[lns2_1_1:_cons])^2
        local reductioncounty = round((1-(`individualcounty' / `emptycounty'))*100,0.1)
        local reductionlocal = round((1-(`individuallocal' / `emptylocal'))*100,0.1)

        display "Reduction county: `reductioncounty' "
        display "Reduction locality: `reductionlocal'"


}



**************************************
* models with contextual vars added  *
**************************************


* code Berlin as eastern 
recode rregion (0=1 East) (1=0 West) (2=1) , gen(macroregion)
drop rregion
ren macroregion rregion
lab var rregion "Macro region"

* Cultural threat 
mixed cultthreat i.male i.higheduc i.agecat i.worker i.unemployedsick  c.localism##c.regexclusion i.kreistyp shareonbenefits shareforeign youngwomen lifeexp i.rregion || district: || location:
est store cultthreat_plus_contextuals 


* Islamophobia
mixed islamophobia i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion i.kreistyp shareonbenefits shareforeign youngwomen lifeexp i.rregion  || district: || location:
est store isphobia_plus_contextuals 


* RWA
mixed rwa i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion i.kreistyp shareonbenefits shareforeign  youngwomen lifeexp i.rregion  || district: || location:
est store rwa_plus_contextuals


* Populism
mixed populism i.male i.higheduc i.agecat i.worker i.unemployedsick c.localism##c.regexclusion i.kreistyp shareonbenefits shareforeign  youngwomen lifeexp i.rregion  || district: || location:
est store populism_plus_contextuals

* Regional exclusion
mixed regexclusion i.male i.higheduc i.agecat i.worker i.unemployedsick i.kreistyp shareonbenefits shareforeign  youngwomen lifeexp i.rregion || district: || location:
est store regex_plus_contextuals

mixed localism i.male i.higheduc i.agecat i.worker i.unemployedsick i.kreistyp shareonbenefits shareforeign  youngwomen lifeexp i.rregion || district: || location:
est store localism_plus_contextuals


* Make a table 

esttab   cultthreat_plus_contextuals isphobia_plus_contextuals rwa_plus_contextuals populism_plus_contextuals conv_plus_contextuals  using multi-level-models.tex, drop(0.male 0.higheduc 1.agecat 0.worker 0.unemployedsick 1.kreistyp 0.rregion) label transform(ln*: exp(2*@) 2*exp(2*@))  varlabels(lns1_1_1:_cons "Variance: county" lns2_1_1:_cons "Variance: locality" lnsig_e:_cons "Variance: person" _cons Constant) eqlabels(none)  cells(b(star fmt(3))) booktabs alignment(D{.}{.}{-1})  replace

* Also save as rtf/Word
esttab   cultthreat_plus_contextuals isphobia_plus_contextuals rwa_plus_contextuals populism_plus_contextuals conv_plus_contextuals  using multi-level-models.rtf, drop(0.male 0.higheduc 1.agecat 0.worker 0.unemployedsick 1.kreistyp 0.rregion) label transform(ln*: exp(2*@) 2*exp(2*@))  varlabels(lns1_1_1:_cons "Variance: county" lns2_1_1:_cons "Variance: locality" lnsig_e:_cons "Variance: person" _cons Constant) eqlabels(none)  cells(b(star fmt(3)))  replace


* Looks like the best way to insert additional midrules is ... sed
shell sed -i '16s/$/\\\midrule/' multi-level-models.tex
shell sed -i '26s/$/\\\midrule/' multi-level-models.tex

* Do marginsplots for the five models
* The values for regexclusion correspond to p10 p25 p50 p75 p90

set scheme plotplain 

est restore cultthreat_plus_contextuals 
margins, at(regexclusion=(1.6/5.6) localism=(3.5 6 7) ) 
marginsplot, subtitle("Cultural threat")
graph export margins-full-cultthreat.pdf, replace 

est restore isphobia_plus_contextuals 
margins, at(regexclusion=(1.6/5.6) localism=(3.5 6 7) ) 
marginsplot, subtitle("Islamophobia")
graph export margins-full-islamophobia.pdf, replace

est restore populism_plus_contextuals 
margins, at(regexclusion=(1.6/5.6) localism=(3.5 6 7) ) 
marginsplot, subtitle("Populism")
graph export margins-full-populism.pdf, replace

est restore rwa_plus_contextuals 
margins, at(regexclusion=(1.6/5.6) localism=(3.5 6 7) ) 
marginsplot, subtitle("Submission/Aggression")
graph export margins-full-rwa.pdf, replace


* RTF Table for regexclusion as dependent variable 
esttab   regex_plus_contextuals  using regexclusion-multi-level-model.rtf, drop(0.male 0.higheduc 1.agecat 0.worker 0.unemployedsick 1.kreistyp 0.rregion) label transform(ln*: exp(2*@) 2*exp(2*@))  varlabels(lns1_1_1:_cons "Variance: county" lns2_1_1:_cons "Variance: locality" lnsig_e:_cons "Variance: person" _cons Constant) eqlabels(none)  cells(b(star fmt(3)))  replace




* I think this saves only the last model :-(
est save multilevel-models, replace 

* OK. Now create a new dataset containing BLUPs for the four variables at the district level
est restore cultthreat_plus_contextuals
predict cultthreatblup, reffects relevel(district)

est restore isphobia_plus_contextuals
predict islamophobiablup, reffects relevel(district)

est restore rwa_plus_contextuals
predict rwablup, reffects relevel(district)

est restore populism_plus_contextuals
predict populismblup, reffects relevel(district)

est restore conv_plus_contextuals
predict convblup, reffects relevel(district)

est restore regex_plus_contextuals
predict regexblup, reffects relevel(district)




preserve
collapse *blup,by(district)
rename district GKZ
save district-blups, replace 
restore 
